Unjamming dynamics: the micromechanics of a seismic fault model 
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The unjamming transition of granular systems is investigated in a seismic fault model via three 
dimensional Molecular Dynamics simulations. A two-time force-force correlation function, and a 
susceptibility related to the system response to pressure changes, allow to characterize the stick-slip 
dynamics, consisting in large slips and microslips leading to creep motion. The correlation function 
unveils the micromechanical changes occurring both during microslips and slips. The susceptibility 
encodes the magnitude of the incoming microslip. 
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In a number of industrial processes and natural phe- 
nomena, such as earthquakes or landslides, disordered 
solid granular systems start to flow. This solid-to-liquid 
transition, known as unjamming, occurs either on de- 
creasing the confining pressure P, or increasing the ap- 
plied shear stress cr. Understanding the properties of 
this transition is a big challenge due to the absence of 
an established theoretical framework for granular mate- 
rials. A proposed analogy with the glass transition [ij of 
thermal systems has recently triggered the study of the 
jamming transition via numerical investigations of sys- 
tems of soft frictionless particles at zero applied shear 
stress where the only control parameter is the pres- 
sure (or the density). As the unjamming transition is 
approached by decreasing the confining pressure, the vi- 
brational spectrum develops an excess of low frequency 
modes, known as soft-modes, leading to the identifica- 
tion of a length scale which diverges on unjamming f^. 
This length scale is related to the emergence of an in- 
creasingly heterogeneous response as the system moves 
towards the transition A different approach to the 
study of the unjamming transition has been followed in 
a two dimensional numerical study [4] and in a number 
of experiments [B-S], where the applied shear stress is 
controlled via a spring mechanism, as the one in Fig. [T^. 
A stick-slij3 motion characterized by a complex slip size 
statistics ^7] is recovered at high confining pressures P 
and small driving velocities Vd- This stick-slip dynamics 
is altered by the presence of noise 0. Analogous results 
have been found at fixed strain rate 13, 

In this Letter we tackle this problem via three di- 
mensional Molecular Dynamics simulations of a model 
of a seismic fault (Fig . [T K ), where grains play the role 
of the gouge [l-d, [ifl El- Numerical details are given 
in [13, • The micromecanical mechanisms leading to 
the transition are analyzed at a level of spatial and tem- 
poral resolution not considered before. 

Stick-slip dynamics - For the investigated values 
of the parameters [13], the system is characterized by 
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FIG. 1: (Color online) (a) The system consists of granular par- 
ticles (light grains) confined between two rigid plates (dark 
grains) at constant pressure. The top plate is driven via a 
spring mechanism, where an extremum of the spring is at- 
tached to the plate and the other is pulled at constant veloc- 
ity. See the supplementary materials available online for an 
animation of a slip event in the real and in the force space [T^ . 
(b) Time evolution of the top plate position X{t). The inset 
indicates that between two large slips there are many mi- 
croslips. (c) Configuration of normal force network. Stronger 
forces correspond to thicker and darker lines, (d) Slip size 
distribution for slips and microslips. 



slip begins and ends when the velocity of the top plate 
becomes, respectively, larger and smaller than a small 
enough threshold. We measure the displacements AA" of 
the top plate due to slips, and compute their distribution 
n{/S.X) (Fig. [TJi). For slips smaller than ~ O.IL^,, where 
Lx is the system length, the distribution follows a power 
law, n(AA) oc AA~^ with /3 ~ 1.85, in agreement with 
experimental values for earthquakes [l6j . Larger slips are 
almost periodic in time and roughly follow a lognormal 
distribution with a characteristic size A A ~ Q.QLx- Sum- 
marizing, the dynamics consists in the occurrence of al- 
most periodic large events, here called slips, and of creep 
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FIG. 2: (Color online) The top plate position in a slip event 
X{t). The vertical dashed line indicates the unjamming time 
tu. Inset: Top plate position for two system replicas, which 
evolve with Vd ~ 0. The replica made at time tu remains in 
the jammed configuration (continuous line) , whereas the one 
made at time tu + St slips (dashed line). 



motion characterized by smaller events, here called mi- 
croslips, in agreement with previous experiments 0, [lH . 

Onset of a slip - To understand the mechanisms acting 
at the onset of a slip, we need to the identify its precise 
starting time i„. Here we describe the analysis performed 
on the particular slip occurring at time t ~ 1600. We 
consider a replica of the system at time t, and follow its 
time evolution at zero driving velocity Vd = 0. If the 
replica made at time t resists to the applied stress, then 
t < tu- Conversely, if a slip is observed, t > t^- We 
define tu as the largest time where no slip occurs, and 
we identify it (one for each slip) with an accuracy equal 
to the time step of integration of the equation of motion 
6t [l3| (Fig- HI ■ This procedure is equivalent to a quasi- 
static simulation [l7| around the un-jamming time t„ and 
gives the value of the shear stress above which the system 
starts to flow. 

We have performed a number of checks which suggest 
that no structural changes occur at t„. For instance, the 
comparison of the state of the system at time t„ with 
the one at shortly earlier and later times, shows that no 
contact breaks at tu- We have also considered the distri- 
bution of the parameter A = |ft|//j,|f„| where ft and f„ 
are the tangential and the normal forces. When A > 1 a 
contact breaks as the Coulomb condition is violated. The 
maximum of the probability distribution P{X) gradually 
moves toward 1 as tu is approached, indicating the weak- 
ening of the solid However, neither at t„ the number 
of contacts with A ~ 1 overcomes a given fraction, nor 
they appear to be spatially organized. The absence of 
structural changes at i„ supports a scenario in which the 
system is located in an energy minimum which slowly 
becomes an inflection point at time t„, and therefore the 



smallest eigenvalue of the dynamical matrix continuously 
decreases to zero [H, [l^ • 

Evolution in the force space - When the system sticks 
and the shear stress increases, no macroscopic motion is 
observed. However, the system microscopically changes 
since it sustains an increasing shear stress. The time evo- 
lution of the top plate position X{t) (Fig. [3^), consists in 
an elastic deformation, where X{t) increases very slowly 
in time, interrupted by sudden microslips. To charac- 
terize the evolution of the system in the force space, we 
introduce a two time force-force correlation function for 
the normal forces, defined as 
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where the sum running over all couples of particles {i,j) 
corresponds to a spatial average. An equivalent defini- 
tion holds for the correlation of tangential forces, C*^. 
Being interested in the unjamming transition of the slip 
event shown in Fig. [2l here we fix to — tu, and con- 
sider the evolution of the correlation function C"(i„,t) 
for earlier times, t < tu (Fig. [^la). The force correla- 
tion function C" (and C'^, not shown) exhibits small 
jumps in correspondence of microslips (Fig. [3^), reveal- 
ing the unusual occurrence of bursts in the reorganiza- 
tion of the force network. During these bursts, the en- 
ergy due to the tangential interaction decreases, whereas 
the one due to the normal interaction increases. A pos- 
sible interpretation is in terms of a two force network 
scenario, in which the applied stress a is supported by 
a stress (t„ due to the normal force network, and by a 
stress at due to the tangential forces, a = an + crt- In 
a burst, few contacts break, leading to a decrease of at, 
at ^ at — 5a- A microscopic slips is observed since the 
normal forces quickly adapt and succeed in sustaining the 
applied stress, an ^ an + 5a- This scenario is supported 
by the inset of Fig. [3b, which shows both C" and C*^ 
across a microslip, with t^ its starting time. C" slightly 
increases and overcomes 1, while C*^ exhibits a sharp 
drop due to the breaking of several contacts. 

The force correlation function (Eq. [1]) gives also in- 
sights into the system evolution during a slip. In Fig.|4]we 
plot C"'*^(t„, t) for t > tu, and for comparison the scaled 
top plate position {X{t)—X{tu))/{X{too)—X{tu)), where 
too is a time following the slip event, whose precise value 
does not influence our results. We first notice that the 
forces evolve on a timescale much shorter than the plate 
motion. For instance the correlation functions reach the 
value 0.1, denoting an almost complete relaxation, when 
the top plate moved only by 10% of its total displace- 
ment. The presence of different time scales in the relax- 
ation process is evidenced by the self-scattering correla- 



tion function F(q, t) 



E,-exp[iq- (rj(i) -rj(tu))] 



where Vj represents the position of the jth particle. Since 
the system is sheared along x and confined along z, we 
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FIG. 3: (Color online) Time evolution of the top plate posi- 
tion (a) and of the correlation function of the normal forces 
C"{t,tu) (b). Dashed lines identify the unjamming time tu- 
The inset show C" and C*® across a microslip at t = 1579.5. 



have considered wave vectors along y, q = (0,g,0). For 
large q, F{q,t) probes small scale relaxation, and coin- 
cides with C"{tu,t), as shown in Fig. Conversely, at 
small q, F{q, t) relaxes on a time scale comparable to 
that of the upper plate motion. The relaxation time r^, 
F{q,Tq) — 1/e, indeed increases as q decreases. More- 
over, tangential forces decorrelate before normal ones. 
This can be explained considering the unjamming tran- 
sition as a buckling-like instability of the chains of large 
normal forces, which are sustained by weaker tangential 
contacts. When the weaker sustaining contacts break, 
either the normal forces adapt to sustain the extra load, 
leading to a microslips, or a buckling-like instability oc- 
curs, giving rise to a slip. The same quantities can 
be used to investigate the subsequent jamming transi- 
tion. To this end, the force network at time t is com- 
pared with the force network after the slip event study- 
ing C"'*^(i, too)- Tangential forces correlate after normal 
ones, and makes the force network stable. 

Response to perturbations: slips versus microslips - 
The correlation length ^ in equilibrium systems is mea- 
sured from the response to an external perturbation. The 
susceptibility, for instance, scales like ^^-jj near critical 
points, where r/ is the correlation function critical expo- 
nent. Here we measure the response of the system to a 
pressure change at the onset of both slips and microslips. 
More precisely, at each time t we stop the external drive, 
setting Vd = 0, and introduce a perturbation in the exter- 
nal pressure P, fixed to P' = P(l — a) for a time interval 
Stperb = 0.1. The response Xa at time t is defined as 



Xa{t) = lim 
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where {rf} and {r,^} are the asymptotic states of the per- 



FIG. 4: (Color online) Correlation functions at the unjam- 
ming (upper panels) and at the subsequent jamming transi- 
tion (lower panels) . The left panels show the normal and tan- 
gential force correlation functions C"'*^(fo,t) (symbols) dur- 
ing a slip event. The dashed line shows the time evolution 
of the plate position scaled between and 1. The plain lines 
show the self-scattering function F{q, t) for q — 3, 5, 9, 17. 
The corresponding are shown in the right panels, to is fixed 
equal to tu for the unjamming transition, and to to ~ 1614 
for the subsequent jamming. 



turbed of the unperturbed systems. In the unjammed 
phase, the susceptibility is divergent. In the jammed 
phase, it measures the size of the region of correlated 
particles that respond to the external perturbation, pro- 
viding an estimate of the correlation length. It is a static 
quantity since the time t only indicates the instant at 
which the perturbation is applied. Xa can be also de- 
fined without setting = 0, provided that the charac- 
teristic response time of the system is much smaller than 
the timescale over which the applied stress varies. For a 
wide range of a, the response of the system at times far 
from tu is linear in the perturbation, since Xa does not de- 
pend on a (Fig. [5]). In particular, Xa gradually increases 
in time and drops in correspondence to microslips. The 
inset shows that the microslip size AX depends on Xa 
evaluated just before the slip, as AX ~ Xa- This indi- 
cates that the size of a microslip is already encoded in 
the system state. In fact, considering that Xa increases 
on approaching a microslip, the measured Xa provides a 
lower bound for the magnitude of the incoming event. 

As the unjamming time is approached, the response 
is no longer linear. Xa remains roughly constant until 
it abruptly increases at a time which depends on a, the 
sooner the greater a. This increase is consistent with a 
power-law divergence. Accordingly at each time, namely 
at each value of the applied shear stress, there is a min- 
imum value of the perturbation intensity for slip trig- 
gering. This behavior is in line with the existence of 
a minimum threshold amplitude in the deformation as- 
sociated with seismic waves for earthquake remote trig- 
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FIG. 5: (Color online) Time evolution of the susceptibility 
Xa for different a. (Inset) Size of a microslip AX versus the 
corresponding value Xa- The straight line is Xa, with b ~ 1.2. 



gering j20l]. A difference between slips and microslips is 
in how different particles contribute to Xa in tlie sum 
of Eq. [21 Indeed, tliese contributions are very similar 
for microslips, and liiglily heterogeneous for slips. The 
presence of an heterogeneous response is consistent with 
previous numerical results found at cr = [Sj. 

In conclusion, the absence of precise structural changes 
at the unjamming time, and the bursts observed in 
the prior dynamics, suggest that the increasing exter- 
nal stress progressively modifies the underlying energy 
landscape. 

At each time the system is in an equilibrium position 
which can be seen as local energy minimum of an effec- 
tive energy landscape which depends on the applied shear 
stress. Microslips occur when the local energy minimum 
flattens down as the applied shear stress increases, letting 
the system fall in a neighbor minimum. If there are no 
close minima, a slip occurs, and the system jumps to a 
far away configuration. The deforming energy landscape 
picture also suggests that soft-modes could be found not 
only when unjamming is approached at cr = by de- 
creasing the volume fraction Q, but also [ilj when 
unjamming is approached increasing a. 
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